Rain in Australia
Adrian Romero Hernandez, Enrique Macip Belmonte, Maria Ruiz Teixidor
19 de diciembre de 2019
Introducción y definición de objetivos
Link repositorio en Github Australian Rain
Los datos obtenidos de Kaggle tienen la información atmosférica de varios años en Australia, datos en los que viene incluida desde la velocidad del viento hasta la temperatura, en total 23 variables recogidas durante varios años en distintas ciudades de Australia con las que se obtienen un dataset de 140.000 líneas aproximadamente.
Nuestro objetivo será predecir la lluvia del día siguiente a partir los datos metereológicos del día (vamos a realizar una regresión logística). Nos interesará tener una matriz de confusión dónde los verdaderos positivos y los verdaderos negativos sean altos. Sobretodo nos interesará que los verdaderos positivos, es decir predecir que llueve el día siguiente, lo predice bien el modelo.
Además, no sólo nos interesa predecir bien si llueve o no el siguiente día, sino entender bien qué factores metereológicos tienen más influencia en la respuesta final. Cómo afecta cada uno de ellos según si su valor aumenta o disminuye.
Para ello, antes que nada hay que entender bien las variables del dataset, qué describen, en qué rangos se mueven, sus valores atípicos y cómo se podrían transformar para poder introducirlas al modelo con la distribución óptima.
Con este dataset tan grande lo primero que nos planteamos fue centrarnos en dos cosas:
Utilizar una zona concreta de Australia, sacada de la variable localización, de la cual elegimos 4 ciudades situadas en la costa sureste de Australia.
Utilizar la variable temporal de alguna forma, ya que considerabamos que tenía importancia pero no podíamos usar cada día del año como un dato diferente, por lo que decidimos obtener a partir de la fecha la estación del año en la que estaba cada línea.
## Dimensiones dataset train: 9824 25## Dimensiones dataset test: 1228 25## Dimensiones dataset validación: 1228 25Antes que nada, visualizamos información básica de las ciudades elegidas y estaciones y cómo están relacionadas con la lluvia.
Realizamos un conteo del número de días que han llovido o no en cada una.
Si filtramos por los días en los que ha llovido hoy (RainToday = 1), vemos las veces que ha llovido el día siguiente o no.
Vemos que en Adelaide, Canberra y Melbourne no suele ser fuecuente que llueva el día siguiente si ha llovido hoy. En Sydney en cambio, es más fequente que llueva si ha llovido el día anterior.
Representamos ahora los días que han llovido en función de las ciudades y las estaciones del año.
En Adelaide y Melbourne las estaciones sí influyen más en la frecuencia de días que llueven, mientras que en Sydney y Canberra suele ser más homogénero.
EDA variables individuales
Resumen de las variables.
## Date Season Location
## Min. :2007-11-01 Length:9824 Length:9824
## 1st Qu.:2010-05-14 Class :character Class :character
## Median :2012-08-31 Mode :character Mode :character
## Mean :2012-10-06
## 3rd Qu.:2015-01-23
## Max. :2017-06-25
##
## MinTemp MaxTemp Rainfall Evaporation
## Min. :-8.00 Min. : 4.10 Min. : 0.000 Min. : 0.000
## 1st Qu.: 7.90 1st Qu.:17.10 1st Qu.: 0.000 1st Qu.: 2.400
## Median :11.60 Median :21.40 Median : 0.000 Median : 4.200
## Mean :11.44 Mean :22.01 Mean : 2.165 Mean : 5.018
## 3rd Qu.:15.50 3rd Qu.:26.10 3rd Qu.: 0.600 3rd Qu.: 6.800
## Max. :33.90 Max. :45.80 Max. :119.400 Max. :43.400
## NA's :9 NA's :8 NA's :181 NA's :2493
## Sunshine WindGustDir WindGustSpeed WindDir9am
## Min. : 0.00 Length:9824 Min. : 11.00 Length:9824
## 1st Qu.: 4.10 Class :character 1st Qu.: 31.00 Class :character
## Median : 7.90 Mode :character Median : 39.00 Mode :character
## Mean : 7.15 Mean : 40.64
## 3rd Qu.:10.20 3rd Qu.: 48.00
## Max. :13.90 Max. :106.00
## NA's :2667 NA's :1104
## WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## Length:9824 Min. : 0.00 Min. : 0.00 Min. : 5.00
## Class :character 1st Qu.: 7.00 1st Qu.:13.00 1st Qu.: 56.00
## Mode :character Median :11.00 Median :19.00 Median : 68.00
## Mean :13.27 Mean :18.82 Mean : 67.22
## 3rd Qu.:19.00 3rd Qu.:24.00 3rd Qu.: 80.00
## Max. :63.00 Max. :65.00 Max. :100.00
## NA's :207 NA's :199 NA's :71
## Humidity3pm Pressure9am Pressure3pm Cloud9am
## Min. : 4.00 Min. : 982.3 Min. : 985.5 Min. :0.000
## 1st Qu.:37.00 1st Qu.:1013.6 1st Qu.:1011.4 1st Qu.:2.000
## Median :48.00 Median :1018.6 Median :1016.3 Median :6.000
## Mean :48.95 Mean :1018.5 Mean :1016.2 Mean :4.694
## 3rd Qu.:60.00 3rd Qu.:1023.5 3rd Qu.:1021.1 3rd Qu.:7.000
## Max. :99.00 Max. :1040.2 Max. :1037.8 Max. :9.000
## NA's :30 NA's :201 NA's :196 NA's :4111
## Cloud3pm Temp9am Temp3pm RainToday
## Min. :0.000 Min. :-1.30 Min. : 3.70 Length:9824
## 1st Qu.:2.000 1st Qu.:11.70 1st Qu.:15.90 Class :character
## Median :5.000 Median :15.40 Median :20.00 Mode :character
## Mean :4.696 Mean :15.55 Mean :20.56
## 3rd Qu.:7.000 3rd Qu.:19.30 3rd Qu.:24.50
## Max. :8.000 Max. :38.60 Max. :44.70
## NA's :4298 NA's :24 NA's :17
## RISK_MM RainTomorrow
## Min. : 0.000 Length:9824
## 1st Qu.: 0.000 Class :character
## Median : 0.000 Mode :character
## Mean : 2.164
## 3rd Qu.: 0.800
## Max. :119.400
## Analizamos las variables individuales por separado con distintos gráficos.
En este dataset hay muchos pares de variables que están fuertemente relacionadas, por ejemplo la temperatura máxima y mínima de un día, o la presión a las 9 de la mañana y la presión a las 3 de la tarde. Por ello, en el análisis individual de variables se estudiarán a la vez por una mejor comprensión.
MinTemp y MaxTemp
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La temperatura mínima se podría asimilar a una distribución normal. La temperatura máxima tiene una cola a su derecha en la que aparecen mútliples valores atípicos.
Temp9am y Temp3pm
Las variables Temp9am y Temp3pm son muy parecidas a las temperaturas mínimas y máximas respectivamente.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Estas variables tienen un comportamiento parecido a la temperatura mínima y máxima respectivamente.
Pressure9am y Pressure3pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Parecen tener una distribución normal ambas variables.
Humidity9am y Humidity3pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Cloud9am y Cloud3pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de una variable categórica que indica la fracción del cielo que está cubierta por nubes a las 9am. Está medida en “oktas”, que son unidades de octavos. Una medida cero indica que el cielo está completamente despejado y un ocho indica que está completamente cubierto.
La variable a las diferentes horas se distribuye de forma muy parecida.
WindSpeed9am y WindSpeed3pm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La velocidad del viente suele tomar valores mayores por la tarde, a las 3pm comparado con la mañana.
WindDir9am y WindDir3pm
Las direcciones en las que el viento se mueve son muy diferentes a las 9 y a las 3. Predominan la dirección Norte y la Oeste por la mañana, mientras que por la tarde suele ser más homogéneo.
Velocidad más fuerte del viento
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de la dirección y la velocidad del viento más fuerte en las últimas 24 horas. Predomina el viento en la dirección norte con una media de veolcidad de 40km/h.
Analizamos cuatro variables que no están, a priori, relacionadas por pares.
Rainfall
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Es la cantidad de lluvia que se ha registrado en el día en mm. Como la mayoría de días no llueve, aparecen muchos valores nulos en el histograma que dificulta la normalización de la variable.
Evaporation
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Se trata de la evaporación en mm en las últimas 24h a las 9am. Contiene un gran número de valroes atípicos.
Horas de sol
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Sorprende el número de días en el que las horas de sol es 0, es decir, está nublado todo el día o con rayos de sol muy débiles.
EDA multivariante
Analizamos la relación unas variables con otras.
Gracias a este gráfico demostramos que los pares de variables citados sí están fuertemente relacionados. Por ejemplo, la presión a las 9 de la mañana con la presión a las 3 de la tarde: si una aumenta, la otra también. Destacar que también hay correlaciones inversas: cuando aumenta la variable Sunshine, disminuye la Cloud9am.
Analizamos en general la relación entre las variables con las estaciones, la variable de salida (RainTomorrow) y las ciudades. Para ello seleccionamos una de las variables de los pares (ya que contienen información similar) y el resto. La variable Rainfall no la mostramos pues su distribución es difícil de visualizar. Se analizarán posteriormente con las transformaciones.
Relaciones por estaciones.
La Temperatura (MaxTemp, Temp9am), la presión (Pressure3pm) y evaporación muestran un claro comportamiento diferente según la estación.
Relaciones por ciudades.
La temperatura es la variable dónde se puede observar más claramente que tiene un comportamiento diferente para cada ciudad.
Relaciones por RainTomorrow, si llueve o no.
Las variables interesantes a analizar con más detalle son: la temperatura, la humedad, la presión y los rayos de sol.
A partir de estas relaciones, indagamos con más detalle las relaciones que parecen interesantes.
Comprobamos que, efectivamente, las cuatro variables de temperatura son muy parecidas, como puede observarse en sus distribuciones:
Al estar relacionadas y tener un comportamiento similar posteriormente se estudiará introdudir al modelo interacciones entre éstas.
Veamos cómo se comporta una de ellas según las estaciones y ciudades:
La Temperatura claramente está condicionada por las estaciones y por las ciudades. No hay mucha diferencia entre otoño y primavera. Dada esta información se podría plantear realizar un modelo para cada estación o para cada ciudad.
Comprobamos que los pares de temperaturas están fuertemente relacionados y que por lo tanto será interesante estudiar interacciones entre ellas.
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Analizamos otras variables más en detalle. Veamos cómo se comporta la presión según las estaciones y ciudades:
La presión se diferencia bastante con cada estación y también se aprecia que está bastante relacionada con si llueve o no: la presión suele ser menor cuando llueve.
Comprobamos que las variables de presión están fuertemente relacionadas entre ellas. 
Veamos cómo se comporta la humedad según las estaciones y ciudades:
Vemos que la humedad contiene valores diferentes según si llueve o no y se diferencia básicamente entre otoño e invierno. Efectiavmente la humedad es mayor los días que llueve.
Comprobamos cómo de relacionadas están las variables de humedad entre ellas.
Veamos cómo se comporta la evaporación según las estaciones y ciudades:
Se observan comportamientos diferentes de la evaporación sobretodo en función de las estaciones.
Veamos cómo se comportan las horas de sol según las estaciones y ciudades:
Claramente la variable sunshine se diferencia más en función de la variable de salida RainTomorrow.
Veamos cómo se comporta la velocidad del viento según las estaciones y ciudades:
La velocidad del viento suele ser mayor los días que llueve.
Veamos cómo se comportan las nubes del viento según las estaciones y ciudades: 
Las variables de nubes a las diferentes horas se comportan de manera muy parecida.
Veamos cómo se comportan la velocidad del viento y su dirección según las estaciones y ciudades:
No parece haber mucha diferencia entre la dirección del viento y la velocidad en función de si llueve o no.
Añadir estos gráficos (Rainfall) después de analizar las transformaciones.
Interacción entre variables y combinaciones
Se puede observar en la gráfica siguiente que ninguna variable individual es mejor que otra para diferenciar si llueve o no llueve, estudiaremos MinTemp, MaxTemp, Temp9am, Temp3pm, Pressure9am, Pressure3pm, Humidity9am, Humidity3pm.
Ahora crearemos unas nuevas varibles combinándolas entre sí. Las usaremos en el modelo más adelante para ajustarlo.
Imputación de datos faltantes
Vamos a realizar un estudio de los datos faltantes para después imputarlos utilizando las técnicas junto al resto de parámetros.
Para analizar los posibles datos faltantes y las relaciones utilizamos dos bibliotecas, visdat y naniar. Con la base de datos de train vemos como se reparten los datos faltantes.
## Date Season Location MinTemp MaxTemp
## 0 0 0 9 8
## Rainfall Evaporation Sunshine WindGustDir WindGustSpeed
## 181 2493 2667 1105 1104
## WindDir9am WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## 709 229 207 199 71
## Humidity3pm Pressure9am Pressure3pm Cloud9am Cloud3pm
## 30 201 196 4111 4298
## Temp9am Temp3pm RainToday RISK_MM RainTomorrow
## 24 17 181 0 0
## AvgTemp Temp Pressure Humidity WindSpeed
## 4 2 178 9 180Con estos datos podemos deducir que las variables con más datos faltantes serían (de más a menos datos faltantes): - Cloud3pm - Cloud9am - Sunshine - Evaporation - WindGustDir - WindGustSpeed - WindDir9am
Con estos datos vamos a estudiar cuáles son las variables con más relaciones de los valores faltantes.
Con esta gráfica podemos observar que los datos faltantes más relaccionados son sobretodo Cloud3pm y Cloud9am, esto puede deberse a que estas variables son categóricas y no pueden medirse directamente, ya que necesitan de la evaluación del cielo en ese día y esto puede que este automatizado.
A estas variables se le suman las variables Evaporation y Sunshine, que también están relaccionadas tanto entre sí como con las mencionadas anteriormente.
Por último recalcar que la variable WindGustDir tiene muchos datos faltantes en los que no tiene relación con las antes mencionadas. Esto puede deberse a que miden cosas completamente diferentes y no necesariamente dependen de las otras.
Con todos estos datos y el EDA con el que hemos empezado podemos comenzar a imputar datos faltantes:
Las primeras variables que vamos a evaluar para completar sus datos faltantes son aquellas variables numéricas que se miden en dos momentos del día, esto nos puede ayudar para completar una con la otra ya que en gran medida estas variables son bastante dependientes y podemos sacar una buena correlación para completar una con la otra.
Para esto utilizaremos el metodo de imputación KNN, con estos datos consideramos que se pueden completar perfectamente estos valores faltantes.
Variables Temp9am & Temp3pm:
#Variables Temp9am & Temp3pm: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Temp9am, Temp3pm) %>% marginplot()
Viendo esta gráfica podemos decir que tiene sentido utilizar la imputación mediante KNN.
#Posible imputación simple mediante un algoritmo de clustering k-NN (con la función VIM:kNN)
#Parece que tiene sentido:
train %>% select(Temp9am, Temp3pm) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")
train_imputed <- kNN(train, variable=c("Temp9am","Temp3pm"))
#test_imputed <- kNN(test, variable=c("Temp9am","Temp3pm"))Comprobamos si nuestra muestra inicial modifica su densidad para asegurarnos de que esta imputación funciona y no hemos generado ningún valor que no coincida.
#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Temp9am, na.rm = T), col=2, main="Temp9am")
lines(density(train_imputed$Temp9am), col=3)
plot(density(train$Temp3pm, na.rm = T), col=2, main="Temp3pm")
lines(density(train_imputed$Temp3pm), col=3) 
Como podemos comprobar, apenas se aprecia la linea roja que esta debajo de la densidad inicial por lo que podemos decir que no cambia.
Variables MinTemp & MaxTemp:
#Variables MinTemp & MaxTemp: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(MinTemp, MaxTemp) %>% marginplot()
#Parece que tiene sentido:
train %>% select(MinTemp, MaxTemp) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")
train_imputed <- kNN(train, variable=c("MaxTemp","MinTemp"))
#test_imputed <- kNN(test, variable=c("MaxTemp","MinTemp"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$MaxTemp,na.rm = T),col=2,main="MaxTemp")
lines(density(train_imputed$MaxTemp),col=3)
plot(density(train$MinTemp,na.rm = T),col=2,main="MinTemp")
lines(density(train_imputed$MinTemp),col=3)
En este caso la imputación de datos faltantes también funciona.
Variables Pressure9am & Pressure3pm:
#Variables Pressure3pm & Pressure9am: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Pressure3pm, Pressure9am) %>% marginplot()
#Parece que tiene sentido:
train %>% select(Pressure3pm, Pressure9am) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")
train_imputed <- kNN(train_imputed, variable=c("Pressure3pm","Pressure9am"))
#test_imputed <- kNN(test_imputed, variable=c("Pressure3pm","Pressure9am"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Pressure3pm,na.rm = T),col=2,main="Pressure3pm")
lines(density(train_imputed$Pressure3pm),col=3)
plot(density(train$Pressure9am,na.rm = T),col=2,main="Pressure9am")
lines(density(train_imputed$Pressure9am),col=3)
Podemos utilizar esta imputación en la presión.
Variables Humidity9am & Humidity3pm:
#Variables Humidity9am & Humidity3pm: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Humidity3pm, Humidity9am) %>% marginplot()
#Parece que tiene sentido:
train %>% select(Humidity3pm, Humidity9am) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")
train_imputed <- kNN(train_imputed, variable=c("Humidity3pm","Humidity9am"))
#test_imputed <- kNN(test_imputed, variable=c("Humidity3pm","Humidity9am"))#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Humidity3pm,na.rm = T),col=2,main="Humidity3pm")
lines(density(train_imputed$Humidity3pm),col=3)
plot(density(train$Humidity9am,na.rm = T),col=2,main="Humidity9am")
lines(density(train_imputed$Humidity9am),col=3)
Podemos utilizar esta imputación en la humedad.
Variables Rainfall & Raintoday: Podríamos pensar que estas dos variables van de la mano, ya que si la premisa se cumple, cuando la variable Rainfall es 0 es que ese mismo día no ha llovido, por ello haremos un estudio de los datos faltantes.
Esto quiere decir que los datos faltantes son exactamente los mismos, por lo que vamos a proceder a completar los datos faltantes de la variable rainfall y con ellos completaremos los datos de Raintoday.
Otras variables que puede tener sentido relacionar es la cantidad de lluvia con la humedad.
Variables Humidity9am & Rainfall:
#Variables Humidity9am & Rainfall: método: imputación simple: KNN
par(mfrow=c(1,1))
train %>% select(Humidity9am, Rainfall) %>% marginplot()
#Parece que tiene sentido:
train %>% select(Humidity9am, Rainfall) %>% VIM::kNN() %>% marginplot(., delimiter="_imp")
train_imputed <- kNN(train_imputed, variable=c("Humidity9am","Rainfall"))
# test_imputed <- kNN(test_imputed, variable=c("Humidity9am","Rainfall"))train_imputed[c("Temp9am_imp","Temp3pm_imp" ,"MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Rainfall_imp", "Humidity3pm_imp", "Pressure3pm_imp","Pressure9am_imp")] <- list(NULL)
imputacion_knn_test<-function(train_imputed, test){
train_imputed["dataset"] <- 'train'
test["dataset"] <- 'test'
dataset_knn = rbind(train_imputed, test)
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Temp9am","Temp3pm"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("MaxTemp","MinTemp"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Humidity3pm","Humidity9am"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Pressure3pm","Pressure9am"))
dataset_knn_imputed <- kNN(dataset_knn, variable=c("Humidity9am","Rainfall"))
train_imputed = subset(dataset_knn_imputed, dataset == 'train')
test_imputed = subset(dataset_knn_imputed, dataset == 'test')
test_imputed[c("Temp9am_imp","Temp3pm_imp" ,"MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Rainfall_imp", "Humidity3pm_imp", "Pressure3pm_imp","Pressure9am_imp","dataset")] <- list(NULL)
return (test_imputed) }
test_imputed = imputacion_knn_test(train_imputed, test)#Comprobaciones:
par(mfrow=c(1,2))
plot(density(train$Humidity9am,na.rm = T),col=2,main="Humidity9am")
lines(density(train_imputed$Humidity9am),col=3)
plot(density(train$Rainfall,na.rm = T),col=2,main="Rainfall")
lines(density(train_imputed$Rainfall),col=3)
## Date Season Location MinTemp MaxTemp
## 0 0 0 0 0
## Rainfall Evaporation Sunshine WindGustDir WindGustSpeed
## 0 2493 2667 1105 1104
## WindDir9am WindDir3pm WindSpeed9am WindSpeed3pm Humidity9am
## 709 229 207 199 0
## Humidity3pm Pressure9am Pressure3pm Cloud9am Cloud3pm
## 0 0 0 4111 4298
## Temp9am Temp3pm RainToday RISK_MM RainTomorrow
## 24 17 181 0 0
## AvgTemp Temp Pressure Humidity WindSpeed
## 4 2 178 9 180Ahora que ya hemos logrado imputar la variable de Rainfall la utilizaremos para conseguir los datos faltantes de raintoday:
train_imputed$RainToday[which(is.na(train_imputed$RainToday))] <- ifelse(train_imputed$Rainfall == 0, "No", "Yes")
test_imputed$RainToday[which(is.na(test_imputed$RainToday))] <- ifelse(test_imputed$Rainfall == 0, "No", "Yes")
train_imputed$RainToday <- ifelse(train_imputed$RainToday == "Yes", 1, 0)
train_imputed$RainToday = as.numeric(train_imputed$RainToday)
train_imputed$RainTomorrow <- ifelse(train_imputed$RainTomorrow == "Yes", 1, 0)
train_imputed$RainTomorrow = as.numeric(train_imputed$RainTomorrow )
test_imputed$RainToday <- ifelse(test_imputed$RainToday == "Yes", 1, 0)
test_imputed$RainToday = as.numeric(test_imputed$RainToday)
test_imputed$RainTomorrow <- ifelse(test_imputed$RainTomorrow == "Yes", 1, 0)
test_imputed$RainTomorrow = as.numeric(test_imputed$RainTomorrow )Variables Cloud3pm & Cloud9am: Para estas variable sutilizamos la librería mice con el metodo de la imputación múltiple y comprobamos que no varía la densidad de los valores.
par(mfrow=c(1,1))
train_imputed0 <- mice(train_imputed[,c('Cloud3pm','Cloud9am')], seed=2018, print = F, m = 30)
train_imputed1 <- mice::complete(train_imputed0)
test_imputed0 <- mice(test_imputed[,c('Cloud3pm','Cloud9am')], seed=2018, print = F, m = 30)
test_imputed1 <- mice::complete(test_imputed0)
xyplot(train_imputed0, Cloud3pm ~Cloud9am)
par(mfrow=c(1,2))
plot(density(train$Cloud3pm,na.rm = T),col=2,main="Cloud3pm")
lines(density(train_imputed1$Cloud3pm),col=3)
plot(density(train$Cloud9am,na.rm = T),col=2,main="Cloud9am")
lines(density(train_imputed1$Cloud9am),col=3)
train_imputed$Cloud3pm <- train_imputed1$Cloud3pm
train_imputed$Cloud9am <- train_imputed1$Cloud9am
test_imputed$Cloud3pm <- test_imputed1$Cloud3pm
test_imputed$Cloud9am <- test_imputed1$Cloud9amVariable Evaporation
En esta varible utilizaremos la media de cada estación y con esto sustituiremos los valores faltantes.
summer <- filter(train, Season == "summer")
evaporationMeanSummer <- mean(summer$Evaporation,na.rm = TRUE)
winter <- filter(train, Season == "winter")
evaporationMeanWinter <- mean(winter$Evaporation,na.rm = TRUE)
spring <- filter(train, Season == "spring")
evaporationMeanSpring <- mean(spring$Evaporation,na.rm = TRUE)
fall <- filter(train, Season == "fall")
evaporationMeanFall <- mean(fall$Evaporation,na.rm = TRUE)
evaporationMeanSummer## [1] 7.992145## [1] 5.715113## [1] 4.140863## [1] 2.447918Aquí podemos observar como la variable de evaporación varia bastante en función de la época del año, en este caso por orden tenemos Verano, Primavera, Otoño e Invierno de mayor a menor, como en el resto de datos de esta variable sólo utiliza un decimal haremos lo mismo para nuestras imputaciones.
evaporationMeanSummer <- round(evaporationMeanSummer,1)
evaporationMeanSpring <- round(evaporationMeanSpring,1)
evaporationMeanFall <- round(evaporationMeanFall,1)
evaporationMeanWinter <- round(evaporationMeanWinter,1)
evaporationMeanSummer## [1] 8## [1] 5.7## [1] 4.1## [1] 2.4train_imputed$Evaporation[which(is.na(train_imputed$Evaporation))] <- ifelse(train_imputed$Season == "summer", evaporationMeanSummer,
ifelse(train_imputed$Season == "winter", evaporationMeanWinter,
ifelse(train_imputed$Season == "spring", evaporationMeanSpring,
ifelse(train_imputed$Season == "fall", evaporationMeanFall,
0))))
test_imputed$Evaporation[which(is.na(test_imputed$Evaporation))] <- ifelse(test_imputed$Season == "summer", evaporationMeanSummer,
ifelse(test_imputed$Season == "winter", evaporationMeanWinter,
ifelse(test_imputed$Season == "spring", evaporationMeanSpring,
ifelse(test_imputed$Season == "fall", evaporationMeanFall,
0))))Variable Sunshine
En esta variable también utilizaremos la media de cada estación y con esto sustituiremos los valores faltantes.
SunshineMeanSummer <- mean(summer$Sunshine,na.rm = TRUE)
SunshineMeanWinter <- mean(winter$Sunshine,na.rm = TRUE)
SunshineMeanSpring <- mean(spring$Sunshine,na.rm = TRUE)
SunshineMeanFall <- mean(fall$Sunshine,na.rm = TRUE)
SunshineMeanSummer## [1] 8.630411## [1] 7.671015## [1] 6.725528## [1] 5.67726Aquí podemos observar como la variable de Sunshine varía bastante en función de la época del año, en este caso por orden tenemos Verano, Primavera, Otoño e Invierno de mayor a menor, como en el resto de datos de esta variable solo utiliza un decimal haremos lo mismo para nuestras imputaciones.
SunshineMeanSummer <- round(SunshineMeanSummer,1)
SunshineMeanSpring <- round(SunshineMeanSpring,1)
SunshineMeanFall <- round(SunshineMeanFall,1)
SunshineMeanWinter <- round(SunshineMeanWinter,1)
SunshineMeanSummer## [1] 8.6## [1] 7.7## [1] 6.7## [1] 5.7train_imputed$Sunshine[which(is.na(train_imputed$Sunshine))] <- ifelse(train_imputed$Season == "summer", SunshineMeanSummer,
ifelse(train_imputed$Season == "winter", SunshineMeanWinter,
ifelse(train_imputed$Season == "spring", SunshineMeanSpring,
ifelse(train_imputed$Season == "fall", SunshineMeanFall,
0))))
test_imputed$Sunshine[which(is.na(test_imputed$Sunshine))] <- ifelse(test_imputed$Season == "summer", SunshineMeanSummer,
ifelse(test_imputed$Season == "winter", SunshineMeanWinter,
ifelse(test_imputed$Season == "spring", SunshineMeanSpring,
ifelse(test_imputed$Season == "fall", SunshineMeanFall,
0))))Variables WindDir & WindSpeed: Para imputar los datos faltantes en la dirección del viento y en la velocidad utilizaremos la velocidad media que hemos obtenido antes en función de la dirección que tenía el aire a esa misma hora, así nos aseguramos de dar un valor mas acorde a la realidad con los datos de los que disponemos.
Para ello, primero tenemos que hacer un estudio de las variables y sus datos faltantes.
trainSelect <- train %>% select(WindDir9am, WindDir3pm, WindSpeed9am,WindSpeed3pm)
gg_miss_upset(trainSelect)
Con esta información podemos concluir que los datos que mas faltan son los de la mañana y destacando los datos de la dirección a las 9 de la mañana. Teniendo en cuenta que en el caso de WindDir3pm, WindSpeed9am, WindSpeed3pm, la mayoría de datos NA corresponden con valores faltantes en los demas valores no seremos capaces de imputar dichos valores con estos datos.
## # A tibble: 17 x 2
## WindDir9am avg_WindSpeed9am
## <chr> <dbl>
## 1 E 9.95
## 2 ENE 9.96
## 3 ESE 10.6
## 4 N 19.4
## 5 NE 9.89
## 6 NNE 10.5
## 7 NNW 14.6
## 8 NW 15.8
## 9 S 13.3
## 10 SE 11.2
## 11 SSE 13.4
## 12 SSW 12.7
## 13 SW 13.2
## 14 W 15.9
## 15 WNW 14.0
## 16 WSW 14.8
## 17 <NA> NAPara poder completar todos los valores vamos a imputar datos en la velocidad del viento. Para ello haremos uso del mice con el que gracias a los valores de la presion a la misma hora imputaremos los valores faltantes.
par(mfrow=c(1,1))
train_imputed0 <- mice(train_imputed[,c('Pressure9am','WindSpeed9am')], seed=2018, print = F, m = 30)
train_imputed1 <- mice::complete(train_imputed0)
test_imputed0 <- mice(test_imputed[,c('Pressure9am','WindSpeed9am')], seed=2018, print = F, m = 30)
test_imputed1 <- mice::complete(test_imputed0)
xyplot(train_imputed0, Pressure9am ~WindSpeed9am)
par(mfrow=c(1,1))
train_imputed2 <- mice(train_imputed[,c('Pressure3pm','WindSpeed3pm')], seed=2018, print = F, m = 30)
train_imputed3 <- mice::complete(train_imputed2)
test_imputed2 <- mice(test_imputed[,c('Pressure3pm','WindSpeed3pm')], seed=2018, print = F, m = 30)
test_imputed3 <- mice::complete(test_imputed2)
xyplot(train_imputed2, Pressure3pm ~WindSpeed3pm)
par(mfrow=c(1,2))
plot(density(train$WindSpeed9am,na.rm = T),col=2,main="WindSpeed9am")
lines(density(train_imputed1$WindSpeed9am),col=3)
plot(density(train$WindSpeed3pm,na.rm = T),col=2,main="WindSpeed3pm")
lines(density(train_imputed3$WindSpeed3pm),col=3)
train_imputed$WindSpeed9am <- train_imputed1$WindSpeed9am
train_imputed$WindSpeed3pm <- train_imputed3$WindSpeed3pm
test_imputed$WindSpeed9am <- test_imputed1$WindSpeed9am
test_imputed$WindSpeed3pm <- test_imputed3$WindSpeed3pma <- data.frame(WindSpeed9am= train_imputed$WindSpeed9am, WindSpeed3pm= train_imputed$WindSpeed3pm)
train_imputed$WindGustSpeed[which(is.na(train_imputed$WindGustSpeed))] <- rowMeans(a, na.rm=TRUE)
b <- data.frame(WindSpeed9am= test_imputed$WindSpeed9am, WindSpeed3pm= test_imputed$WindSpeed3pm)
test_imputed$WindGustSpeed[which(is.na(train_imputed$WindGustSpeed))] <- rowMeans(b, na.rm=TRUE)Transformaciones de variables cualitativas
train_imputed[c("Temp9am_imp","Temp3pm_imp" ,"Pressure9am_imp","Pressure3pm_imp","MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Humidity3pm_imp","Rainfall_imp")] <- list(NULL)
test_imputed[c("Temp9am_imp","Temp3pm_imp" ,"Pressure9am_imp","Pressure3pm_imp","MaxTemp_imp","MinTemp_imp", "Humidity9am_imp","Humidity3pm_imp","Rainfall_imp")] <- list(NULL)Transformaciones de variables cuantitativas
La única que no parece seguir una distribución normal es la variable RainFall. Analizamos cual es la mejor transformación que podemos aplicar.
p1 <- train_imputed %>% select(Rainfall, Season) %>%
na.omit() %>%
ggplot(aes(x=Rainfall, colour=Season)) +
geom_density()
aux <- filter(train_imputed,train_imputed$Rainfall>0)
p2 <- aux %>% mutate(inv_Rainfall = 1/ (Rainfall + 1)) %>%
select(Season, inv_Rainfall) %>%
na.omit() %>%
ggplot(aes(x=inv_Rainfall, colour=Season)) +
geom_density()
p3 <- train_imputed %>% mutate(log10_Rainfall = log10(Rainfall + 1)) %>%
select(Season, log10_Rainfall) %>%
na.omit() %>%
ggplot(aes(x=log10_Rainfall, colour=Season)) +
geom_density()
grid.arrange(p1, p2, p3, nrow = 1)
train_imputed$Rainfall_Inv <- 1/ (train_imputed$Rainfall + 1)
test_imputed$Rainfall_Inv <- 1/ (test_imputed$Rainfall + 1)Con esta información podemos concluir que para transformarla lo mejor que podemos hacer es utilizar la inversa y sumarle 1 para poder realizar bien la transformación.
Comprobamos algunas variables que sigan distribución normal.
#Resto de variables: siguen distribución normal
qqnorm(train_imputed$MinTemp, ylab="MinTemp")
qqline(train_imputed$MinTemp, col="red")
Resumen de valores nulos
## [,1]
## Date 0.00000000
## Season 0.00000000
## Location 0.00000000
## MinTemp 0.00000000
## MaxTemp 0.00000000
## Rainfall 0.00000000
## Evaporation 0.00000000
## Sunshine 0.00000000
## WindGustDir 11.24796417
## WindGustSpeed 0.00000000
## WindDir9am 7.21701954
## WindDir3pm 2.33102606
## WindSpeed9am 0.00000000
## WindSpeed3pm 0.00000000
## Humidity9am 0.00000000
## Humidity3pm 0.00000000
## Pressure9am 0.00000000
## Pressure3pm 0.00000000
## Cloud9am 0.00000000
## Cloud3pm 0.00000000
## Temp9am 0.24429967
## Temp3pm 0.17304560
## RainToday 0.00000000
## RISK_MM 0.00000000
## RainTomorrow 0.00000000
## AvgTemp 0.04071661
## Temp 0.02035831
## Pressure 1.81188925
## Humidity 0.09161238
## WindSpeed 1.83224756
## Rainfall_Inv 0.00000000El resto de nulos los eliminamos por ser un pequeño porcentaje.
Modelo 1: selección con Lasso sobre el dataset original.
Estandarización:
Aplicamos la estandarización con las medias y varianzas de train.
#Estandarización
numeric_vars <- c("MinTemp","MaxTemp", "Rainfall", "Evaporation", "Sunshine", "WindGustSpeed",
"WindSpeed9am", "WindSpeed3pm", "Humidity9am", "Humidity3pm","Pressure9am" ,
"Pressure3pm" , "Temp9am", "Temp3pm", "Rainfall_Inv")
preprocessParams <- preProcess(train_imputed[numeric_vars], method=c("center", "scale"))
train_imputed[numeric_vars] <- predict(preprocessParams, train_imputed[numeric_vars])
#Aplicamos la estandarización con las medias y varianzas de train
test_imputed[numeric_vars] <- predict(preprocessParams, test_imputed[numeric_vars])Creación variables dummies para Location, Season, WindGustDir y borrado de las originales.
train_imputed <- dummy_cols(train_imputed, select_columns = c("Location", 'Season', 'WindGustDir', 'WindDir9am', 'WindDir3pm'))
test_imputed <- dummy_cols(test_imputed, select_columns = c("Location", 'Season', 'WindGustDir','WindDir9am', 'WindDir3pm'))
train_imputed$Location <- NULL
train_imputed$Season <- NULL
train_imputed$Date <- NULL
train_imputed$WindGustDir <- NULL
train_imputed$RISK_MM <- NULL
test_imputed$Location <- NULL
test_imputed$Season <- NULL
test_imputed$Date <- NULL
test_imputed$WindGustDir <- NULL
test_imputed$WindDir3pm_NA <- NULL
test_imputed$WindDir9am_NA <- NULL
test_imputed$RISK_MM <- NULLSelección de variables usando Lasso
El modelo utilizado es una regresión logística, realizamos la selección de variables mediante Lasso. Para ello, primero seleccionamos el lambda óptimo y con éste seleccionamos las variables.
train_imputed_sinNA = na.omit(train_imputed)
x_train = train_imputed_sinNA[, !names(train_imputed_sinNA) %in% c("RainTomorrow")]
y_train = train_imputed_sinNA$RainTomorrow
test_imputed_sinNA = na.omit(test_imputed)
x_test = test_imputed_sinNA[, !names(test_imputed_sinNA) %in% c("RainTomorrow")]
y_test = test_imputed_sinNA$RainTomorrow
x = model.matrix(RainTomorrow~ ., train_imputed_sinNA)[,-1]
lambda_seq <- 10^seq(2, -2, by = -.1)
cv.out <- cv.glmnet(x, y_train, alpha = 1, lambda = lambda_seq)
plot(cv.out)
# identifying best lamda
best_lam <- cv.out$lambda.min
lasso_best <- glmnet(x, y_train, alpha = 1, lambda = best_lam)
c<-coef(lasso_best,s=best_lam,exact=TRUE)
inds<-which(c!=0)
variables<-row.names(c)[inds]Variables seleccionados por Lasso
## [1] "(Intercept)" "Rainfall" "Sunshine"
## [4] "WindGustSpeed" "WindDir9amNNE" "WindDir9amS"
## [7] "WindDir3pmNNW" "WindSpeed3pm" "Humidity3pm"
## [10] "Pressure3pm" "Cloud3pm" "RainToday"
## [13] "WindSpeed" "Location_Adelaide" "Season_summer"
## [16] "WindGustDir_NNW" "WindGustDir_ENE" "WindDir9am_NW"
## [19] "WindDir9am_S" "WindDir9am_NNE" "WindDir3pm_NNW"
## [22] "WindDir3pm_NNE"Modelo usando las variables Lasso
new_train = train_imputed_sinNA %>% select("Sunshine","WindGustSpeed", "WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "Humidity3pm", "Pressure3pm", "Cloud3pm", "RainToday", "Location_Adelaide", "Location_Sydney", "Season_summer", "WindGustDir_ENE", "WindGustDir_NNW", "WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "RainTomorrow")
new_test = test_imputed_sinNA %>% select("Sunshine","WindGustSpeed", "WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "Humidity3pm", "Pressure3pm", "Cloud3pm", "RainToday", "Location_Adelaide", "Location_Sydney", "Season_summer", "WindGustDir_ENE", "WindGustDir_NNW","WindDir9am_N", "WindDir9am_NNE", "WindDir9am_S", "WindDir3pm_NNE", "RainTomorrow")
#Modelo
#Entrenamiento
glm_model_train = glm(RainTomorrow~ ., data=new_train, family= binomial)
#Test
glm_test = predict(glm_model_train, newdata = new_test, type = "response")Resumen del modelo
##
## Call:
## glm(formula = RainTomorrow ~ ., family = binomial, data = new_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6951 -0.5643 -0.3426 -0.1585 2.8729
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.39044 0.10181 -23.479 < 2e-16 ***
## Sunshine -0.50184 0.03863 -12.990 < 2e-16 ***
## WindGustSpeed 0.47714 0.04123 11.573 < 2e-16 ***
## WindDir9am_N 0.09919 0.09587 1.035 0.300826
## WindDir9am_NNE 0.39545 0.13169 3.003 0.002674 **
## WindDir9am_S -0.33364 0.14500 -2.301 0.021398 *
## WindDir3pm_NNE 0.49834 0.17700 2.816 0.004869 **
## Humidity3pm 0.91951 0.04239 21.690 < 2e-16 ***
## Pressure3pm -0.51749 0.03935 -13.152 < 2e-16 ***
## Cloud3pm 0.05059 0.01401 3.610 0.000306 ***
## RainToday 0.55232 0.07221 7.649 2.03e-14 ***
## Location_Adelaide 0.63438 0.08212 7.725 1.12e-14 ***
## Location_Sydney 0.22857 0.08438 2.709 0.006754 **
## Season_summer -0.16317 0.08425 -1.937 0.052778 .
## WindGustDir_ENE -0.45050 0.18304 -2.461 0.013847 *
## WindGustDir_NNW 0.34854 0.12695 2.745 0.006043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8823.1 on 8168 degrees of freedom
## Residual deviance: 6117.7 on 8153 degrees of freedom
## AIC: 6149.7
##
## Number of Fisher Scoring iterations: 5Vemos que las variables más significativas son los rayos de sol, la velocidad del viento, la humedad, la presión, si ha llovido hoy y si es la ciudad Adelaide, dónde llueve mucho más respecto el resto de ciudades en invierno.
Curva ROC del modelo
Matriz de confusión
#Matriz de confusión
umbral_dec = 0.46
glm_test <- ifelse(glm_test >= umbral_dec, 1, 0)
glm_test <- factor(glm_test, levels = c(0,1))
tabla_conf <- table(glm_test, new_test$RainTomorrow)
caret::confusionMatrix(tabla_conf, positive = '1')## Confusion Matrix and Statistics
##
##
## glm_test 0 1
## 0 731 115
## 1 59 111
##
## Accuracy : 0.8287
## 95% CI : (0.8041, 0.8514)
## No Information Rate : 0.7776
## P-Value [Acc > NIR] : 3.216e-05
##
## Kappa : 0.4569
##
## Mcnemar's Test P-Value : 3.052e-05
##
## Sensitivity : 0.4912
## Specificity : 0.9253
## Pos Pred Value : 0.6529
## Neg Pred Value : 0.8641
## Prevalence : 0.2224
## Detection Rate : 0.1093
## Detection Prevalence : 0.1673
## Balanced Accuracy : 0.7082
##
## 'Positive' Class : 1
## Observando la matriz de confusión se puede ver que el modelo no obtiene buenos resultados pues el número de falsos positivos es muy alto y casi igual al de verdaderos positivos. Por lo que el modelo predice bien cuando no llueve pero no acaba de predecir bien cuando sí va a llover.
Probamos otros modelos para ver si el resultado mejora.
Modelo 2: selección con Lasso sobre el dataset con interacción de variables.
Como hay muchas variables que están muy relacionadas entre ellas probamos a combinarlas para introducirlas al modelo y que de esta manera no se pierda información.
Para este nuevo modelo usaremos las nuevas variables AvgTemp, Temp, Pressure, Humidity creadas de las medias de las variables a las 9am y a las 3pm.
#Estandarización
numeric_vars <- c("MinTemp","MaxTemp", "Rainfall", "Evaporation", "Sunshine", "WindGustSpeed",
"WindSpeed9am", "WindSpeed3pm", "Humidity9am", "Humidity3pm","Pressure9am",
"Pressure3pm", "Temp9am", "Temp3pm", "Rainfall_Inv", "AvgTemp", "Temp", "Pressure", "Humidity", "WindSpeed" )
preprocessParams <- preProcess(train_imputed[numeric_vars], method=c("center", "scale"))
train_imputed[numeric_vars] <- predict(preprocessParams, train_imputed[numeric_vars])
#Aplicamos la estandarización con las medias y varianzas de train
test_imputed[numeric_vars] <- predict(preprocessParams, test_imputed[numeric_vars])train_imputed$WindSpeed9am <- NULL
train_imputed$WindSpeed3pm <- NULL
train_imputed$Humidity9am <- NULL
train_imputed$Humidity3pm <- NULL
train_imputed$MinTemp <- NULL
train_imputed$MaxTemp <- NULL
train_imputed$Pressure9am <- NULL
train_imputed$Pressure3pm <- NULL
train_imputed$Temp9am <- NULL
train_imputed$Temp3pm <- NULL
test_imputed$WindSpeed9am <- NULL
test_imputed$WindSpeed3pm <- NULL
test_imputed$Humidity9am <- NULL
test_imputed$Humidity3pm <- NULL
test_imputed$MinTemp <- NULL
test_imputed$MaxTemp <- NULL
test_imputed$Pressure9am <- NULL
test_imputed$Pressure3pm <- NULL
test_imputed$Temp9am <- NULL
test_imputed$Temp3pm <- NULLSelección de variables usando Lasso
train_imputed_sinNA = na.omit(train_imputed)
x_train = train_imputed_sinNA[, !names(train_imputed_sinNA) %in% c("RainTomorrow")]
y_train = train_imputed_sinNA$RainTomorrow
test_imputed_sinNA = na.omit(test_imputed)
x_test = test_imputed_sinNA[, !names(test_imputed_sinNA) %in% c("RainTomorrow")]
y_test = test_imputed_sinNA$RainTomorrow
dim(test_imputed_sinNA)## [1] 1020 73## [1] 8186 75#Primero vemos Lasso
x = model.matrix(RainTomorrow~ ., train_imputed_sinNA)[,-1]
lambda_seq <- 10^seq(2, -2, by = -.1)
cv.out <- cv.glmnet(x, y_train, alpha = 1, lambda = lambda_seq)
plot(cv.out)
# identifying best lamda
best_lam <- cv.out$lambda.min
lasso_best <- glmnet(x, y_train, alpha = 1, lambda = best_lam)
c<-coef(lasso_best,s=best_lam,exact=TRUE)
inds<-which(c!=0)
variables<-row.names(c)[inds]
variables## [1] "(Intercept)" "Rainfall" "Sunshine"
## [4] "WindGustSpeed" "WindDir9amNE" "WindDir9amNNE"
## [7] "WindDir9amNW" "WindDir9amS" "WindDir9amSSW"
## [10] "WindDir3pmNNE" "WindDir3pmNNW" "WindDir3pmSSW"
## [13] "WindDir3pmWNW" "Cloud3pm" "RainToday"
## [16] "Pressure" "Humidity" "WindSpeed"
## [19] "Location_Adelaide" "Location_Canberra" "Season_spring"
## [22] "WindGustDir_ESE" "WindGustDir_NNW" "WindDir9am_SSW"
## [25] "WindDir9am_NW" "WindDir3pm_WNW" "WindDir3pm_SSW"
## [28] "WindDir3pm_NNW" "WindDir3pm_NNE"new_train = train_imputed_sinNA %>% select("AvgTemp", "Temp", "Pressure", "Humidity", "WindSpeed", "RainTomorrow")
new_test = test_imputed_sinNA %>% select("AvgTemp", "Temp", "Pressure", "Humidity", "WindSpeed", "RainTomorrow")
#Modelo
#Entrenamiento
glm_model_train = glm(RainTomorrow~ ., data=new_train, family= binomial)
#Test
glm_test = predict(glm_model_train, newdata = new_test, type = "response")
#Evaluación modelo
#Summary
summary(glm_model_train)##
## Call:
## glm(formula = RainTomorrow ~ ., family = binomial, data = new_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5023 -0.6633 -0.4295 -0.1834 2.9040
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.56426 0.03436 -45.519 < 2e-16 ***
## AvgTemp 0.31869 0.15059 2.116 0.0343 *
## Temp -0.27327 0.16029 -1.705 0.0882 .
## Pressure -0.69921 0.03584 -19.507 < 2e-16 ***
## Humidity 1.12150 0.04517 24.830 < 2e-16 ***
## WindSpeed 0.19759 0.03291 6.003 1.93e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8839.2 on 8185 degrees of freedom
## Residual deviance: 6993.8 on 8180 degrees of freedom
## AIC: 7005.8
##
## Number of Fisher Scoring iterations: 5#Tabla de ganancia
logistic_gains_table <- blr_gains_table(glm_model_train, data = new_train)
#Curva ROC
blr_roc_curve(logistic_gains_table)
#Matriz de confusión
umbral_dec = 0.46
glm_test <- ifelse(glm_test >= umbral_dec, 1, 0)
glm_test <- factor(glm_test, levels = c(0,1))
tabla_conf <- table(glm_test, new_test$RainTomorrow)
caret::confusionMatrix(tabla_conf, positive = '1')## Confusion Matrix and Statistics
##
##
## glm_test 0 1
## 0 739 142
## 1 54 85
##
## Accuracy : 0.8078
## 95% CI : (0.7823, 0.8316)
## No Information Rate : 0.7775
## P-Value [Acc > NIR] : 0.009957
##
## Kappa : 0.3555
##
## Mcnemar's Test P-Value : 5.156e-10
##
## Sensitivity : 0.37445
## Specificity : 0.93190
## Pos Pred Value : 0.61151
## Neg Pred Value : 0.83882
## Prevalence : 0.22255
## Detection Rate : 0.08333
## Detection Prevalence : 0.13627
## Balanced Accuracy : 0.65318
##
## 'Positive' Class : 1
## Vemos que utilizando la combinación de las variables no mejora el modelo. El número de falsos positivos es mayor que el de verdaderos positivos. El modelo continua prediciendo bien cuando no va a llover pero su predicción cuando llueve no es exacta.
Los próximos pasos serían realizar un modelo para cada una de las ciudades estudiadas o para cada estación.